/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
** ****************************************************************** */

//
// 文件: GPUElementKernels.h
// 创建时间: 2025年10月2日
// 说明: GPU单元刚度矩阵计算核函数声明
// 作者: OpenSees GPU加速项目组
//
// 功能:
//   - 声明GPU单元刚度矩阵计算核函数
//   - 支持Truss、Beam、Quad等常用单元类型
//   - 提供高性能并行刚度矩阵计算
//

#ifndef GPUElementKernels_h
#define GPUElementKernels_h

#include "GPUElementSupport.h"
#include <cuda_runtime.h>

// Forward declaration
struct ConcreteState;

//==============================================================================
// GPU单元刚度矩阵计算核函数声明
//==============================================================================

/**
 * @brief GPU并行计算Truss单元刚度矩阵
 *
 * @param d_elementData 单元数据数组（GPU）
 * @param d_stiffness 单元刚度矩阵输出缓冲区（GPU）
 * @param d_nodeDisplacements 节点位移向量（用于几何非线性，可选）
 * @param numElements 单元数量
 *
 * @note 每个线程处理一个Truss单元
 * @note 刚度矩阵大小：6x6（2节点×3DOF）
 */
__global__ void computeTrussStiffness(
    GPUElementData* d_elementData,
    double* d_stiffness,
    double* d_nodeDisplacements,
    int numElements
);

/**
 * @brief GPU并行计算四边形单元刚度矩阵
 *
 * @param d_elementData 单元数据数组（GPU）
 * @param d_stiffness 单元刚度矩阵输出缓冲区（GPU）
 * @param d_nodeDisplacements 节点位移向量（用于几何非线性，可选）
 * @param numElements 单元数量
 *
 * @note 每个线程处理一个Quad单元
 * @note 刚度矩阵大小：8x8（4节点×2DOF）
 * @note 使用2x2高斯积分
 * @note 支持平面应力/平面应变
 */
__global__ void computeQuadStiffness(
    GPUElementData* d_elementData,
    double* d_stiffness,
    double* d_nodeDisplacements,
    int numElements
);

/**
 * @brief GPU并行计算Beam单元刚度矩阵
 *
 * @param d_elementData 单元数据数组（GPU）
 * @param d_stiffness 单元刚度矩阵输出缓冲区（GPU）
 * @param d_nodeDisplacements 节点位移向量（用于几何非线性，可选）
 * @param numElements 单元数量
 *
 * @note 每个线程处理一个Beam单元
 * @note 2D梁刚度矩阵大小：6x6（2节点×3DOF）
 * @note 3D梁刚度矩阵大小：12x12（2节点×6DOF）
 */
__global__ void computeBeamStiffness(
    GPUElementData* d_elementData,
    double* d_stiffness,
    double* d_nodeDisplacements,
    int numElements
);

//==============================================================================
// GPU辅助设备函数声明
//==============================================================================

/**
 * @brief 计算2D四边形单元的Jacobian矩阵
 *
 * @param coords 单元节点坐标[4][2]
 * @param xi 自然坐标系xi值
 * @param eta 自然坐标系eta值
 * @param J Jacobian矩阵输出[2][2]
 * @param detJ Jacobian行列式输出
 *
 * @return 0-成功，-1-失败（奇异矩阵）
 */
__device__ int computeQuadJacobian(
    double coords[4][2],
    double xi,
    double eta,
    double J[2][2],
    double& detJ
);

/**
 * @brief 计算2D四边形单元的形函数导数（物理坐标系）
 *
 * @param xi 自然坐标系xi值
 * @param eta 自然坐标系eta值
 * @param invJ Jacobian逆矩阵[2][2]
 * @param dN_dx 形函数对x的导数输出[4]
 * @param dN_dy 形函数对y的导数输出[4]
 */
__device__ void computeQuadShapeFunctionDerivatives(
    double xi,
    double eta,
    double invJ[2][2],
    double dN_dx[4],
    double dN_dy[4]
);

/**
 * @brief 计算平面应力材料矩阵
 *
 * @param E 弹性模量
 * @param nu 泊松比
 * @param D 材料矩阵输出[3][3]
 */
__device__ void computePlaneStressD(
    double E,
    double nu,
    double D[3][3]
);

/**
 * @brief 计算应变-位移矩阵（B矩阵）
 *
 * @param dN_dx 形函数对x的导数[4]
 * @param dN_dy 形函数对y的导数[4]
 * @param B 应变-位移矩阵输出[3][8]
 */
__device__ void computeStrainDisplacementMatrix(
    double dN_dx[4],
    double dN_dy[4],
    double B[3][8]
);

//==============================================================================
// C接口声明（用于从C++代码调用）
//==============================================================================

#ifdef __cplusplus
extern "C" {
#endif

/**
 * @brief 启动Truss单元刚度矩阵计算核函数
 *
 * @param numElements 单元数量
 * @param d_elementData 单元数据（GPU）
 * @param d_stiffness 刚度矩阵缓冲区（GPU）
 * @param d_nodeDisplacements 节点位移（GPU，可选）
 * @param stream CUDA流（用于异步执行）
 *
 * @return 0-成功，-1-失败
 */
int launchTrussStiffnessKernel(
    int numElements,
    const GPUElementData* d_elementData,
    double* d_stiffness,
    const double* d_nodeDisplacements,
    cudaStream_t stream = 0
);

/**
 * @brief 启动Quad单元刚度矩阵计算核函数
 *
 * @param numElements 单元数量
 * @param d_elementData 单元数据（GPU）
 * @param d_stiffness 刚度矩阵缓冲区（GPU）
 * @param d_nodeDisplacements 节点位移（GPU，可选）
 * @param d_materialStates 材料状态数组（GPU）
 * @param d_strainTotal 总应变数组（GPU）
 * @param maxDOFPerElement 单元最大DOF数
 * @param stream CUDA流（用于异步执行）
 *
 * @return 0-成功，-1-失败
 */
int launchQuadStiffnessKernel(
    int numElements,
    const GPUElementData* d_elementData,
    double* d_stiffness,
    const double* d_nodeDisplacements,
    ConcreteState* d_materialStates_committed,
    ConcreteState* d_materialStates_trial,
    double* d_strainTotal,
    double* d_strainCurrent_trial,
    int maxDOFPerElement,
    cudaStream_t stream = 0
);

/**
 * @brief 启动Beam单元刚度矩阵计算核函数
 *
 * @param numElements 单元数量
 * @param d_elementData 单元数据（GPU）
 * @param d_stiffness 刚度矩阵缓冲区（GPU）
 * @param d_nodeDisplacements 节点位移（GPU，可选）
 * @param stream CUDA流（用于异步执行）
 *
 * @return 0-成功，-1-失败
 */
int launchBeamStiffnessKernel(
    int numElements,
    const GPUElementData* d_elementData,
    double* d_stiffness,
    const double* d_nodeDisplacements,
    cudaStream_t stream = 0
);

/**
 * @brief 启动节点坐标更新核函数（Phase 2优化）
 *
 * 计算公式：currentCoord = initialCoord + displacement
 *
 * 优化目标：
 * - 避免CPU端遍历所有节点提取坐标
 * - 在GPU上直接计算当前坐标
 * - 预期收益：消除每步2-3ms CPU开销
 *
 * @param numNodes 节点数量
 * @param d_currentCoords 输出：当前节点坐标（GPU）
 * @param d_initialCoords 输入：初始节点坐标（GPU，首次缓存）
 * @param d_displacements 输入：节点位移（GPU，每步传输）
 * @param stream CUDA流（用于异步执行）
 *
 * @return 0-成功，-1-失败
 */
int launchUpdateNodeCoordinatesKernel(
    int numNodes,
    double* d_currentCoords,
    const double* d_initialCoords,
    const double* d_displacements,
    cudaStream_t stream = 0
);

#ifdef __cplusplus
}
#endif

#endif // GPUElementKernels_h
